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1.  Introduction 

Stochastic  linear  programs  can  be  formulated  for  a  variety  of  applications.  Some  examples  include  airline 
scheduling  (Ferguson  and  Dantzig  [1956]),  financial  planning  (Kusy  and  Ziemba  [1986]),  energy  modeling 
(Birge  [1987])  and  water  resource  planning  (PreTcopa  and  Szantai  [1978]).  The  basic  model  we  consider  here 
is  the  stochastic  linear  program  with  recourse  in  the  following  general  form: 

mmx{cTx  +  Q{x)\Ax  =  b,  x  >  0} 

where 

Q{x)=[Q{x,t,<l>)P{dt,d<l>) 

and  the  recourse  function  is  defined  as 

Q{x,  £,  4>)  =  min{gry|Wy  =  £  -  Tx,  u  +  <f>  >  y  >  0}, 

where  i  6  3Jn\  y  €  3Jn3,  b  6  9Rmi,  and  (£,(f>)  is  a  random  vector  on  the  probability  space  (dlm*+n* ,  7,  P) 
with  support,  Hx$.  The  vectors,  c,  q,  and  u,  and  matrices,  A,  W,  and  T  are  dimensioned  correspondingly. 
The  fundamental  problem  in  stochastic  programming  is  to  evaluate  the  integral  of  Q.  In  this  paper,  we 
describe  a  method  for  finding  an  upper  bound  on  Q  that  requires  a  polynomial  number  of  operations  in  the 
number  of  random  variables. 

Previous  results  in  bounding  expressions  for  Q  are  described  in  Birge  and  Wets  [1986a].  The  bounds  are 
based  on  the  convexity  and  positive  homogeneity  of  Q.  The  first  result  is  due  to  Jensen  [1906] 's  inequality 
which  provides  a  lower  bound  on  Q.  The  usefulness  of  this  lower  bound  is  that  it  requires  an  evaluation  of  Q 
at  one  point  (the  mean  of  the  random  variables)  and  has  been  found  to  be  generally  sharp  in  some  practical 
examples  (see,  e.g.,  Hausch  and  Ziemba  [1983]).  Madansky  [1959]  provided  an  upper  bound  following 
Edmundson  [1956]  that  is  based  on  the  theory  of  moment  spaces  and  amounts  to  weighting  the  extreme 
points  of  the  support  of  the  random  variables.  Ben-Tal  and  Hochman  [1972]  and  Huang,  Ziemba  and  Ben-Tal 
[1977]  refined  this  bound  for  independent  random  variables.  Dupa£ova  [1976]  formulated  a  bound  of  the  same 
general  type  for  dependent  random  variables  that  was  extended  to  unbounded  ranges  and  non-polyhedral 
sets  in  Gassmann  and  Ziemba  [1986].  Frauendorfer  [1986]  provided  a  sharper  bound  in  the  bounded  range, 
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dependent  variable  case,  and  Birge  and  Wallace  [1986]  gave  a  bound  and  method  for  refinement  for  special 
cases  of  dependent  random  variables. 

The  upper  bounds  mentioned  above  all  have  the  property  that  they  are  solutions  to  moment  problems 
with  varying  conditions.  Dupacova's  work  on  minimax  solutions  (2a£kova  [1966])  led  to  these  conclusions 
and  to  the  use  of  the  generalized  moment  problem.  Ermoliev,  et  al.  [1987]  provided  a  general  programming 
framework  for  solving  the  general  problem.  It  is  used  in  Birge  and  Wets  [1987]  for  bounds  with  piecewise 
linear  approximations  on  moment  constraints  and  in  Cipra  [1985]  with  first  and  second  moment  constraints. 

The  problem  with  each  of  these  bounds  is  that  they  require  an  exponentially  increasing  number  of 
function  evaluations  as  the  number  of  random  variables  increases.  An  alternative  for  this  situation  was  given 
by  the  ray  approximation  procedure  in  Birge  and  Wets  [1986a].  This  uses  the  sublinearity  property  of  the 
recourse  function  to  obtain  a  separable  function  that  majorizes  Q.  This  approach  is  generalized  in  Birge 
and  Wets  [1986b].  Wallace  [1987b],  on  the  other  hand,  formulated  a  procedure  that  applies  to  problems  in 
which  the  recourse  function  involves  the  solution  of  a  network  problem.  Our  procedure  is  a  combination  and 
generalization  of  these  two  basic  approaches.  The  algorithm  we  give  provides  a  separable  piecewise  linear 
function  that  bounds  Q  throughout  the  support  of  the  random  variables  and  can  be  easily  evaluated. 

Section  2  presents  our  basic  algorithm  and  the  separable  piecewise  linear  upper  bound  (  SPLU).  Its 
properties  are  described  in  Section  3.  Section  4  gives  an  illustrative  small  example  and  provides  comparison 
with  the  upper  bound  of  Edmundson  and  Madansky.  Extensions  of  the  basic  algorithm  and  conclusions  are 
given  in  Section  5. 

2.  The  Basic  Algorithm 

We  give  a  general  method  for  finding  an  upper  bound  on  the  expected  value  of  the  value  of  a  linear 

program  with  random  right-hand  sides  and  random  upper  bounds  on  the  variables.    To  simplify  notation 

and  to  establish  general  results,  we  consider  the  following  system  : 

Aix  =  &i  +  £ 

A2x  =  b2  (1) 

0  <  x  <  c  +  <j> 
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where  Ax  G  9imiXn,  A2  6  gfc("l-"li)x'l)  (^il^)7  =  ,4  is  the  coefficient  matrix,  (&i|62)T  =  &  is  the  fixed 
part  of  the  right-hand  side,  c  is  the  fixed  part  of  the  bounds  on  the  variables,  £  is  the  random  availability 
of  resources  and  (j>  is  the  random  part  of  the  variable  capacities,  where  <f>  >  0.  We  assume  that  there  is  a 
positive  probability  that  <j>  =  0.  Next  define  Q(£,  (f>)  by 

<?(£,«*)  =  min{gTa:|(l)}.  (2) 

Finally  define  x(£,  <t>,  d~ ,  d+)  as  the  set  of  x-vectors  satisfying 

AyX=£ 

A2x  =  0  (3) 

d~  <  x<  <j>  +  d+. 
Our  goal  is  to  find  an  upper  bound  on  Q(£,<£),  or,  more  precisely,  on  EQ(f,  <f>).    We  do  this  by  finding  a 
separable  piecewise  linear  function  U(£,  <f>)  defined  by 

•=i  urx<-(s-ei)iff<ii 

where  &  =  E&,  and  #((£)  is  a  piecewise  linear  function  in  <f>. 

Algorithm  1 

Step  0:  Find  Q{£,0)  with  optimal  solution  x°,  where 

(efBo  1{b  +  £)     if  t  is  basic, 
0  if  t  is  nonbasic  at  lower  bound, 

c(x)  if  t  is  nonbasic  at  upper  bound. 

Assume  for  simplicity  that  the  first  m  variables  are  basic.     Let  xi+   =   (5,71e,,0,0, . . .  ,0)  and  x*~   = 

(-BQ1ei,0,0,...,Q)  where  t  =  1,2,...  .n^.  Let 

^(i)  =  max-  x°(t)  -  Y,x>+{i)y>+  -^V"**)^- 

y=2  j=a 


subject  to 


yJ+-yJ-      =  £,-£■ 

efn  <^<^max,y  =  2,...,m1. 


Let 


subject  to 


"»l 


mi 


P'ii)  =  min  -  z°(t)  -  Yl  «y+(*V+  -  J2  x°~  (*V~  +  c(z) 

]=2  j=2 


efn  <ei<efax,y=2,...,m1, 

for  all  t  =  1, . . . ,  n. 

If  a1^')  >  0  for  some  i  or  pl{i)  <  0  for  some  i,  let  x*+  =  a1^')  =  (0, ...  ,0),  and  pl{i)  =  x°(i)  +  c(t')  for  all 

i  =  1, . . . ,  mi  and  go  to  Step  1  with  r  =  1. 

Otherwise,  check 

e(i)  =  max  —  x1+(i)y+  —  x1-  (i)y~ 


subject  to 


and 


subject  to 


y+  -  y~  =  £i  -  & 

/•min   ^   c     s'   /-max 

£1     SaS  ?i 


/(»')  =  min  -  i1+(t')y+  -  i1    (t')y 


y+-y~  =$i-  £i 


^rn  <  &  <  ^rax- 


If  a1  +  e  <  0  and  /91  +  /  >  0,  then  Q(£,  <£)  is  linear  in  f ,  go  to  Step  4. 

Otherwise,  let  r  =  1  and  go  to  Step  1. 
Step  1:  If  £™ax  <  +oo,  solve 


min{grx  |  X[(£rmax  -  fr)er,  0,  ar,&]}  =  qTxr+{1?°*  -  £r) 


Else  (let  PZ{i)  =  oo  if  fir(i)  =  +co,  Pl(i)  =  0  otherwise)  and  solve 


T„r  + 


min{qr x  \  X{er,0,0,  PI)}  =  q1  x 


If  £"in  >  -oo,  solve 


min{gTx  |  X[(^  -  Zr)eriO,ar,0r]}  =  qTxr+(-^  +  |r). 


Else  (let  #(t)  =  oo  if  £r(t)  =  +00,  Bl(i)  =  °  otherwise)  and  solve 


T„r- 


mm{qTx  |  x(-«r,  0, 0,  #)}  =  -q1  x 


If  Step  1  was  entered  with  xl±  =  (0, . . . ,  0)  for  all  t,  go  to  Step  2;  otherwise,  go  to  Step  4. 
Step  2:  For  t  =  1 n,  solve 


ar+1(0  =  max-x°(t')-    £    z>+(*V"+  -    £    iJ'~(*V" 


subject  to 


subject  to 


/9-+1(t)  =  min-x°(t)-    X)   *y+(*V+-    £   ^"(tV'-+c(0 


yy+_yy-    =  ^  -  ey 

£fn     <^<€T*,j^r+l. 


Step  S:  If  r  <  m,  let  r  =  r  +  1  and  go  to  Step  1. 
Otherwise,  go  to  Step  4. 
Step  4:  Find 


mi 


mi 


subject  to 


subject  to 


for  i  =  1, . . . ,  n. 


a*(i)  =  max  -  *°(t)  -  J>"+(tV+  "  X>_(*V' 

y=i  y=i 


j^+-yy-    =  &-& 

£fn    <t,<trx,j=l,...,m1, 

mi  n%i 

/?*(*)  =  min  -  x°(t)  -  j2^+(tV+  ~  ]C^'*(*V"  +c(0 

y=i  y=i 


^+-yy-    =  &-& 

efn  <ey<efax,y=i,...,m1, 


Let  x*  =  axgmin  {qT x  |  x(0,^max,  a*,y9*)}.  Find  a  conformal  realization  of  x*  (Rockafellar  [1984, p.  455]),so 
that 

x*  =  YJ  ajkZJt  with  ak  >  0, 

such  that  x*{i)  >  0  =>  x*k(i)  >  0  and  x*(i)  <  0  =>  x£(z)  <  0,  and  x*{i)  =  0  =>  x£(i)  =  0.  An  algorithm  for 
finding  such  a  realization  is  the  "painted  index  algorithm"  in  Rockafellar  [l984,p.476].  Paint  all  columns  Aj 

of  A  such  that 

{white  if  x*  (j)  >  0, 

black  if  x*  (j)  <  0, 

red  if  x*  (j)  =  0. 

Let  k  —  1.  Pivot  until  a  Tucker-tableau  is  reached  in  which  there  is  a  compatible  column.  This  will  always 

be  possible  in  our  case.    Let  the  compatible  column  be  A'-,  and  let  F  be  the  set  of  indices  for  the  basic 

columns  in  the  final  Tucker  tableau.  We  now  have  that 

i€F 

If  Aj  is  white,  let 

(A'^i)    if  ieF, 

x%{%)  =  {  1  ift  =  j, 

v  0  otherwise. 

If  Aj  is  black,  reverse  all  signs  in  x£.  (Note  that  the  sign  convention  in  a  Tucker  tableau  is  opposite  of  the 
convention  in  the  standard  simplex  tableau.) 

Let  afc  =  min{z*(z)/x£(ii),  x£(t)  ^  0},  x*(i)  =  x*(i)  —  akX^{i)  and  re-paint  every  column  for  which  x*(i)  =  0 
red. 

If  x*  7^  0,  let  k  =  k  +  1  and  repeat.  Otherwise,  go  to  Step  5  with  the  conformal  realization  X^fc=i  akxk- 

Step  5:  Using  the  cost  coefficients  qT xx± ,  find  E(U(t;,  <f>).  This  amounts  to  performing  mi  simple  line 
integrals. 

Step  6:  If  x*(i)  >  0  (so  that  z£(t)  >  0,  VA;),  we  are  using  a  variable  x(i)  with  random  capacity  (3*(i)+<f>i{>  4>i). 
If  x*(i)  <  0,  we  are  using  a  variable  x(i)  with  deterministic  capacity  a*(i)(<  0).  We  shall  in  the  following 
assume  that  each  variable  x(z),such  that   x*(t)    ^   0,   has  associated  with  it  a  random  arc  capacity  <f>*. 
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If  x*(i)   <  0,  we  have  Pr{#   =  a*(i)}  =  1,  if  z*(i")  >  0,  <f>*  =  <j>  +  0*{i).    For  each  k  =  l,...,Kt  let 
Ik  =  ]C9(0xfc(0(^  °)-  Sort  the  primal  supports  x*k  such  that  ?i  <  ft  <  ■  ■  •  <  <7jc.  Let  k  =  1,  p  =  0  (where 
p  will  become  EH(<f>)). 
Step  7:  Let  P  =  {*  |  x%{%)  ^  0}.  Consider  the  random  variable 

pk  =  max{O,mmP{<t>;/xl{i)}},0k  €  [0,ofc]. 

Find  Eqicfik-  (This  work  amounts  to  increasing  the  capacity  of  each  conformal  flow  until  the  first  variable 
capacity  is  met.  This  continues  on  each  conformal  flow.  Details  are  given  for  the  network  case  in  Wallace 
[1987b].)  Let  p  =  p  +  Eqk/3k,  and  <f>*  =  <f>*  —  afcz£.  If  k  =  K  or  if  qk+i  =  0,  stop  with  p  =  EH (<£),  otherwise 
let  k  =  k  +  1  and  repeat  Step  7. 

End. 

The  value  obtained  in  Algorithm  1  is  indeed  an  upper  bound  on  the  expected  linear  program  value. 

Theorem.  The  value  SPLU=  E^ft/ff,  (/>)]  obtained  in  Algorithm  1  is  an  upper  bound  on  Q  =  E€i^,[Q(^,  <f>)]. 

Proof:  The  proof  requires  only  showing  that  x  =  x°  +  $I(xJ  +  (£j  —  £)+  +  **"(£  —  ly)  +  )  +  H2{PkQkXk)  is 
feasible  in  x(£,  <£*,0,  c).  This  is  obtained  by  noting  that  the  definitions  of  ir±,  ar,  and  ^9r  in  Steps  0  to  2 
maintain  feasibility  for  <j>*M 

The  algorithm  as  described  above  is  our  basic  version.    We  prove  certain  properties  of  it  in  the  next 
section.  In  Section  5,  we  present  alternative  versions  of  some  of  the  steps  in  Algorithm  1. 

3.  Properties  of  the  Upper  Bound 

The  purpose  of  this  section  is  to  show  that  the  upper  bound  presented  in  this  paper  has  some  desirable 
properties  and  to  relate  the  procedure  to  other  bounding  methods. 

3.1  Exact  Bounds  for  Linear  Problems 

All  other  bounds  used  in  stochastic  programming  are  exact  whenever  Q(£,  <f>)  is  linear  in  f  and  <f>  over 
the  support  of  the  random  variables.  This  is  true  of  the  Madansky  upper  bound,  the  piecewise  linear  upper 
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bound  in  the  pure  network  case  (Wallace  [1987b]),  the  linear  upper  bound  on  the  expected  max  flow  in  a 
network  (Wallace  [1987a]),  the  Jensen  lower  bound,  and  the  sublinear  approximation  in  Birge  and  Wets  (if 
the  random  variables  have  unbounded  support).  All  acceptable  bounds  should  have  this  property. 

Property  1:  The  bound  SPLU  given  by  Algorithm  1  is  exact  if  Q(i,4>)  is  a  linear  function. 

Proof:  Assume  Q(£,  <f>)  is  linear  in  £  and  <f>  and  that  the  reduced  cost  of  a  non-basic  variable  is  always 
different  from  zero  (a  dual  non-  degeneracy  assumption).  Then  E£/(£, <f>)  =  EQ(^,<f>).  Of  course,  if  Q  is 
linear,  it  can  be  written  as 


mi 


Q{i,t)  =  Q&o)  +  Y,  fkitk  -  &)  +  £M;- 

Clearly,  Step  0  provides  us  with  Q(£,  0).  Also,  if  Q  is  linear,  a1  +  e  <  0,  /91  +  /  >  0  in  Step  0, since  the  basis 
corresponding  to  Q(i-,  0)  is  feasible  for  all  f  G  E.  Hence,  fk  =  qxk+  =  —  qxk~ .  Therefore,  if  Q  is  linear,  the 
algorithm  will  discover  the  coefficients  of  £  in  Step  0  and  then  go  to  Step  4. 

Let  us  define  a  variable  i  to  be  stochastic  if  <£™ax  >  0,  otherwise,  it  is  deterministic.  Consider  the 
conformal  realization  of  x*  =  Ylakxk-  First  note  that  z£  is  an  elementary  vector  (Rockafellar  [1984,  p. 453]). 
This  means  that  there  is  no  way  to  split  x£  into  two  or  more  other  vectors  where  at  least  one  has  fewer 
non-zeroes  than  x£. 

Assume  there  exists  an  elementary  vector  y  such  that  y(i)  ^  0  for  more  than  one  stochastic  random 
variable.  Then  fix  the  value  of  <f>i  at  0  for  all  variables  except  for  those  with  y(i)  ^  0.  Then,  Q  would  not  be 
linear.  (Compare  with  the  random  variable  /?  in  Step  7.)  Hence,  if  Q  is  linear,  there  is  no  elementary  vector 
with  more  than  one  stochastic  variable. 

Now,  assume  that  we  have  found  two  elementary  vectors  t/i  and  3/2  ,  such  that  they  share  the  stochastic 
variable  i  (i.e.,  j/i(t)  7^  0,y2(t)  ^  0).  Also  assume  that  q\/yi{i)  ^  ^/sM1)-  (The  variable  ft  defined  as  in 
Step  6.)  Let  all  <j>j  =  0  for  t  ^  j.  Then  Q  is  not  linear  in  variable  i,  because  the  marginal  gain  of  increasing 
<f>i  is  not  the  same  in  both  elementary  vectors.  Hence,  two  elementary  vectors  can  only  share  a  stochastic 
variable  if  <7i/t/i(t)  =  92/j/2(0-  (This  corresponds  to  two  circuits  in  a  pure  network  that  have  the  same  cost 
and  share  an  arc  with  a  random  capacity.)  Of  course,  hi  =  <?i/j/i(0- 

Hence,  if  Q  is  linear,  no  elementary  vector  x£  has  more  than  one  stochastic  variable  and  two  elementary 

9 


vectors  can  only  share  a  stochastic  variable  if  they  have  the  same  cost  (in  the  sense  described  above).  Since 
Step  6  only  creates  elementary  vectors  x£,  the  random  variable  Pk  in  Step  7  is  linear  in  its  single  random 
variable.  Hence,  our  method  produces  the  exact  solution.! 

3.2  The  Bound  is  Polynomial 

The  Edmundson-Madansky  bound  requires  that  Q{£,  <f>)  be  solved  in  all  extreme  cases  of  £  and  <f>.  There 
are  2mi+ni  such  points;  hence,  the  method  is  exponential  in  the  number  of  stochastic  variables.  Only  for 
very  moderate  values  of  ni  and  rr»i  is  it  possible  to  apply  this  bound. 

The  major  goal  of  this  paper  is  therefore  to  find  a  good  upper  bound  that  can  be  computed  in  a  number 
of  operations  that  is  polynomial  rather  than  exponential  in  the  number  of  random  variables. 

Property  2:  Algorithm  1  calculates  SPLU  in  a  number  of  operations  that  is  polynomial  in  the  number  of 
random  variables. 

Proof:  The  amount  of  work  is  in  the  worst  case: 

Step  0:  1  LP  (a1,^1  can  be  found  by  inspection). 

Step  1:  2mi  LPs. 

Step  4:  1  LP  to  find  i*.  The  conformal  realization  is  independent  of  i%i  and  mi.  (The  worst  case  is  n  LP's, 
n  <  ni.) 

Step  5:  The  integration  is  a  constant  amount  of  work  for  each  random  variable. 

Step  7:  Finding  E^/?*  amounts  to  checking  the  mi  *  max{mj}  (in  the  worst  case)  possible  values  of  /?fc.  The 
value  mi  is  the  total  number  of  possible  values  for  fa.  This  has  to  be  done  not  more  than  n  times  (since  the 
number  of  zeroes  increases  by  one  for  each  k). 

Hence,  the  algorithm  is  linear  in  ni  and  rriiM 

3.3  Relation  to  Networks 

The  method  presented  in  this  paper  is  closely  related  to  the  network  method  in  Wallace  [1987b].   The 
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major  difference  is  in  Step  1,  where  we  only  solve  two  networks  in  the  network  case  and  not  2mj  as  here. 

Below  is  a  short  network  interpretation  of  some  of  the  vectors  and  scalars  used  in  the  algorithm  to  help 
in  its  understanding. 
Step  0:  The  variable  i*+  shows  how  the  flow  changes  on  the  basic  arcs  as  the  supply  at  node  i  is  increased 

by  one  unit  or  (the  demand  is  decreased).  Hence, 

{+1     if  arc  j  is  a  forward  arc  on  the  path  from  node  i  to  the  slack  node, 
—  1     if  arc  j  is  a  reverse  arc, 
0        if  arc  j  is  not  on  the  path. 

xl~  is  similarly  defined  for  increased  demand  (or  decreased  supply). 

a1(i)  >  0  implies  that  with  the  chosen  set  of  paths  (x*)  there  are  supply/demand  combinations  that  give  a 
negative  flow  on  arc  t  ,  even  when  we  disregard  node  1. 

/91(i)  >  0  implies  that  with  the  chosen  paths,  there  are  supply/demand  combinations  that  overuse  arc  i  even 
when  we  do  not  consider  node  1. 

Step  1:  xf  are  still  paths,  but  not  along  a  basis.  Both  basic  and  non-basic  arcs  are  used.  If  Step  1  finishes 
successfully,  we  have  actually  replaced  the  original  network  by  a  star-shaped  network  (where  the  slack  node 
is  in  the  center  of  the  star).  The  arc  going  from  the  center  node  to  node  i  has  unit  cost  qx~ ,  the  arc  in 
the  other  direction  has  unit  cost  qxf .  The  way  we  have  used  ar  and  0r  has  guaranteed  that  whatever 
combination  we  get  of  supply  and  demand,  sending  that  flow  along  the  paths  xf-  would  be  feasible  and  cost 
the  same  as  in  the  star-shaped  network. 

Hence,  we  have  found  an  upper  bounding  simple  recourse  problem  (Wets  [1983])  .  In  stochastic  programming 

this  approximation  depends  on  the  actual  value  of  the  first  stage  decisions  (as  in  the  recourse  function  in 

the  introduction).  Hence,  in  some  sense,  it  is  a  local  approximation. 

Step  4:  a*  <  0  shows  how  much  flow  can  be  sent  along  the  original  arcs  in  the  negative  direction  without 

making  that  total  flow  negative  (whatever  the  supply/demand  is).  Similarly,  /?*  shows  how  much  is  left  of 

the  capacity  in  the  arcs  in  the  worst  case. 

x*  is  just  a  circulation  in  the  network,  and  x*k  are  circuits  of  minimal  length  (in  terms  of  the  number  of 

arcs  in  them),  a^  shows  how  much  flow  the  circuit  can  take  (or,  more  precisely,  how  much  flow  it  has  been 

allotted.) 
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Step  7:  /?*  is  a  random  variable  describing  the  capacity  of  circuit  k. 

3.4  Relation  to  Sublinear  Approximations 

The  separable  piecewise  linear  upper  bound  is  also  a  generalization  of  the  ray  function  approximation  in 
Birge  and  Wets  [1986a]  and  its  extension  in  the  sublinear  approximation  in  Birge  and  Wets  [1986b].  These 
procedures  find  the  value  of  Q(£,  <j>)  in  different  coordinate  directions  to  again  obtain  a  separable  function 
that  can  easily  be  integrated.  The  approach  in  Birge  and  Wets  [1986b]  uses  varying  choices  of  the  coordinate 
system  that  leads  to  an  extension  of  the  SPLU  bound  given  here.  This  extension  would  involve  solving  for 
iJ+  and  x3~  in  different  directions  so  that  a  variety  of  bounds  could  be  obtained. 

The  ray  function  approximation  amounts  to  solving  for 

qTx3+  =  min{<7Tx  |  Ax  =  e3,x  >  0} 

and 

q7 ' x?~  =  min{qTx  \  Ax  =  —e3,x  >  0}. 

These  values  of  x3  +  and  x3~  are  then  used  in  U{£,<f>)  as  in  SPLU.  The  extension  is  to  use  the  elements  of 
other  coordinate  systems  in  place  of  ±e3  in  the  definitions  (i.e.,  use  some  vectors  d3  that  form  a  basis  for 
9?n).  This  procedure  can  be  used  in  Algorithm  1  to  obtain  an  alternative  bound. 

The  sublinear  approximation  with  varying  directions  has  been  found  to  produce  accurate  approximations 
in  a  variety  of  examples.  The  advantage  of  the  SPLU  bound  is  that  it  applies  to  bounded  regions  so  it  may 
be  used  on  partitions  of  the  support  of  the  random  variable  in  a  refinement  procedure  in  solving  a  stochastic 
program.  Algorithm  1  also  incorporates  the  procedures  for  handling  random  bounds  that  often  arise  in 
practical  examples. 

3.5  Finiteness 

There  is  no  guarantee  that  our  upper  bound  is  finite,  i.e.,  that  all  linear  programs  that  must  be  solved 
are  feasible.  An  infinite  bound  of  course  results  if  EQ(f ,  <f>)  —  +oo,  i.e.  the  problem  itself  is  infeasible,  but 
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it  can  also  be  that  EQ(£,  <f>)  <  +00,  whereas  EU{£,4>)  =  +00.  This  is  not  always  avoidable.  We  note  that 
the  only  other  polynomial  upper  bound,  the  ray  approximation,  is  never  better  than  our  bound  (assuming 
the  possible  extensions  mentioned  above),  and  that  exponential  bounds  may  be  necessary  in  some  cases. 

3.6  Partitioning 

When  approximations,  such  as  the  one  in  this  paper,  are  used  in  two-stage  stochastic  programming,  a 
comparison  is  made  with  a  lower  bound  (EL)  usually  based  on  Jensen's  inequality.  Then,  if  EC/  -  EL  is 
too  large  (according  to  some  rule),  the  support  rectangle  (for  independent  random  variables)  is  partitioned 
into  smaller  rectangles  called  cells,  and  the  bounding  procedures  are  applied  to  these  cells,  which  in  turn  are 
weighted  by  their  probability. 

Hence,  whenever  a  partition  is  called  for,  one  must  decide  which  cell  to  partition  and,  along  which 
coordinate  direction,  to  perform  the  partition.  With  an  upper  bounding  method  that  can  take  on  the  value 
+00  even  for  a  feasible  problem,  one  should  clearly  partition  the  cell  where  EU  —  +00,  along  the  coordinate 
direction  that  was  being  treated  when  the  infeasibility  was  discovered.  This  provides  a  dynamic  scheme  in 
which  the  algorithm  is  applied  on  each  cell  until  either  an  infinite  value  is  obtained  or  the  difference  between 
lower  and  upper  bounds  is  above  the  acceptable  threshold.  A  partition  is  made  in  either  instance.  Partition 
strategies  are  discussed  in  Birge  and  Wets  [1986a],  Birge  and  Wallace  [1986]  and  Frauendorfer  and  Kail 
[1986]. 

4.  Examples 

In  this  section,  we  first  present  a  small  example  to  illustrate  the  bound.  We  then  give  computational 
results  on  a  larger  problem  from  energy  modeling  (Louveaux[l987]). 

4.1  A  Problem  with  Two  Random  Variables 

The  first  example  is  a  problem  with  two  random  variables  and  without  random  capacities.  We  wish  to 
find  bounds  on  EQ(£)  where 

Q(Z)  =  min  xx  +  x2  +  x3  +  x4  +  10i5  +  10x6  (4.1) 
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subject  to 


xi      +3x2     +13  -*6  =  6     (4-2) 

3ix      +x2  +x4  -xq     =  £2     (4.3) 

X\ , . . . ,  Xq  >  0, 


where  £1  and  £2  are  uniformly  distributed  on  [1,4].  This  problem  is  illustrated  in  Figure  1,  where  A.i  refers 
to  the  ith  column  in  the  constraint  matrix  of  (4.2-4.3).  We  follow  Algorithm  1  step  by  step. 

Step  0: 

(i)  Find  Q(£)  =  1.25. 

i°  =  (0.625,0.625,0,0,0,0). 

xi+  =  (-0.125,0.375, 0,0,0, 0);x1_  =  (0.125,-0.375,0,0,0,0). 

i2+  =  (0.375, -0.125, 0,0, 0,0);  i2_  =  (-0.375,0.125,0,0,0,0). 

ax(l)  =  -0.0625;  a1  (2)  =  -0.4375;(Note  :  £r(t)  =  +00.) 

Now  e(l)  =  0.1875  >  —al(\)  =  0.0625,  so  go  to  Step  1.  Note  in  Figure  1  that  we  have  essentially  moved 
along  the  vertical  line  through  £.  The  bound  a1  recorded  the  (negative)  minimum  multiples  of  the  vectors 
A.i  and  A_2  for  points  along  that  line.  The  value  e(l)  recorded  the  greatest  change  in  the  multiple  of  A,i 
from  the  multiple  for  £  for  other  points  along  the  horizontal  line  through  £.  The  function  is  not  linear 
because  this  change  is  greater  than  the  minimal  multiple  (ai(l))  for  movement  in  the  vertical  direction. 

Step  1. 

Solve  min{9Tz  |  x[1.5ci,0,  (-0.0625,  -0.4375,  0,  0,0,0),  00)]}  =  1.125  =  (0.75)  *  (1.5)  =  qTx1+{^ax  -  £1), 
where  i1+  =  (-0.0625,0.1875, 1.0,0,0,0),  and 

min{qTx  |  x[-l-5ei,  0,  (-0.0625,  -0.4375,0,0, 0, 0),  00)]}  =  1.375  =  (0.9167)  *  (1.5)  =  gTz1-(£i  -  £rin), 

where  i1"  =  (-0.0625,-0.4375,0,0.625,0.125,0).  Next  go  to  Step  4. 

Step  4: 

We  can  skip  this  since  there  are  no  random  bounds.  Step  5  then  is  the  terminal  step. 

Step  5: 
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Here  we  compute  E£/(f)  = 

Qi$  +  f         *1+(£i  "  L)dF(^)  +  f         x1"^  -  ZJdFfa) 
+  /       z2+(6  -  6)^(6)  +  /       *2"(6  -  6)^(6) 

=  1.25  +  0.5(0.75)  (0.75)  +  (0.5)  (0.9167)  (0.75)  +  0.5(0.25)  (0.75)  +  (0.5)(-0.25)(0.75)  =  1.875. 

End. 

So,  we  have  SPLU  =  1.875.  We  compare  this  with  the  Edmundson-Madansky  (EM)  bound.  In  this 
example,  the  EM  bound  assigns  equal  weights  to  the  values  of  Q(£)  at  each  of  the  extreme  points  of  5. 
Hence, 

EM  =  0.25  *  (Q(l,  1)  +  Q(l,  4)  +  Q(4, 1)  +  Q(4,  4))  =  1.625. 

The  EM  bound  is  better  than  the  SPLU  bound  but  this  difference  may  be  eliminated  by  refinements  of  the 
SPLU  bound.  We  describe  possible  refinements  in  Section  5. 

4.2  Computational  Results  for  an  Energy  Model 

The  usefulness  of  the  SPLU  bound  is  best  demonstrated  on  a  practical  example  in  which  the  number  of 
random  variables  varies.  We  wish  specifically  to  observe  the  performance  of  SPLU  relative  to  the  EM  bound 
as  the  number  of  random  variables  increases.  The  performance  is  measured  in  the  sharpness  of  the  bound 
and  the  computational  effort.  As  a  practical  example,  we  consider  the  small  energy  model  in  Louveaux 
[1987].  We  do  not  consider  random  bounds  because  that  is  directly  analogous  to  the  network  case  discussed 
in  Wallace  [1987b]. 

In  this  example,  we  have  four  technologies  which  can  be  used  to  satisfy  three  demands  at  varying  costs. 
High  cost  "backstop"  technologies  are  also  available  to  satisfy  demand  so  the  problem  is  feasible  for  any 
demand  realization.  The  randomness  occurs  in  the  capacity  of  the  technologies  and  the  demands.  This  allows 
from  one  to  seven  random  variables.  The  examples  were  also  chosen  with  varying  ranges  (narrow, medium, 
and  wide)  on  the  random  variables  resulting  in  twenty-one  sets  of  examples.  We  assume  uniform  distributions. 
This  assumption  favors  bounds  (such  as  the  Edmundson-Madansky  bound)  that  place  weights  at  extreme 
values  since  other  distributions  generally  have  more  mass  around  the  center  of  the  support. 
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The  experiments  were  conducted  on  the  Amdahl  5860  at  The  University  of  Michigan  Computing  Center. 
The  SPLU  and  EM  bounds  were  both  implemented  in  FORTRAN  codes  using  the  same  linear  programming 
routine  LPM-1  (Pfefferkorn  and  Tomlin  [1976]).  Each  bound  was  computed  for  each  of  the  twenty-one  test 
problems.  The  Jensen  inequality  lower  bound  was  also  computed  to  determine  the  values  of  the  upper 
bounds  relative  to  the  lower  bounds.  The  results  tire  given  in  Table  1. 

The  results  in  Table  1  show  that  the  polynomial  bound  SPLU  does  not  generally  provide  as  accurate  a 
bound  as  the  EM  bound,  but  that  as  the  number  of  random  variables  increases  the  computational  time  in 
SPLU  increases  much  less  rapidly  than  the  time  for  EM.  In  these  examples,  the  growth  of  time  for  SPLU 
is  indeed  approximately  linear  (gaining  ten  milliseconds  for  each  random  variable),  while  the  time  for  EM 
approximately  doubles  as  each  new  random  variable  is  introduced.  This  demonstrates  that  the  real  advantage 
of  the  SPLU  bound  is  in  problems  with  a  large  number  of  random  variables  where  the  EM  bound  cannot  be 
computed.  These  results  are  comparable  with  the  results  in  Wallace  [1987b]  for  networks. 

Refinements  are  also  possible  to  reduce  the  error  in  SPLU.  In  Section  5,  a  refinement  scheme  using 
parametric  linear  programming  is  introduced.  The  use  of  different  coordinate  directions  is  another  possibility 
as  mentioned  above.  For  unbounded  ranges,  the  resulting  sublinear  approximation  values  were  reduced  up 
to  thirty  percent  from  the  coordinate  direction  values  for  a  similar  set  of  test  problems  (Birge  and  Wets 
[1986b]). 
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PROBLEM1 

TIME2 

VALUE 

EM 

SPLU 

EM 

SPLU 

Jensen 

1-NAR 

9 

10 

182.75 

182.75 

182.75 

1-MED 

10 

18 

220.50 

220.25 

220.00 

1-WID 

10 

21 

385.50 

341.50 

297.50 

2-NAR 

16 

20 

183.38 

183.06 

182.75 

2-MED 

17 

26 

220.50 

220.50 

220.00 

2-WID 

22 

35 

389.85 

389.10 

297.50 

3-NAR 

22 

28 

183.38 

183.38 

182.75 

3-MED 

25 

38 

221.38 

222.50 

220.00 

3-WID 

33 

40 

433.60 

439.58 

297.50 

4-NAR 

51 

41 

184.09 

185.50 

182.75 

4-MED 

44 

41 

227.22 

255.18 

220.00 

4-WID 

47 

45 

434.26 

469.50 

297.50 

5-NAR 

94 

52 

184.19 

186.44 

182.75 

5-MED 

87 

51 

227.41 

278.38 

220.00 

5-WID 

75 

54 

434.35 

499.18 

297.50 

6-NAR 

163 

54 

185.58 

192.49 

182.75 

6-MED 

193 

61 

235.91 

303.35 

220.00 

6-WID 

149 

72 

443.52 

524.30 

297.50 

7-NAR 

329 

64 

186.23 

215.48 

182.75 

7-MED 

366 

70 

236.27 

328.35 

220.00 

7-WID 

304 

76 

444.12 

549.30 

297.50 

1  -  Number  of  random  variables  -  range  of  random  variables. 

2  -  CPU  milliseconds. 

Table  1.  Results  for  Edmundson-Madansky  and  SPLU  bounds. 
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5.  Extensions  and  Conclusions 

The  SPLU  bound  can  be  refined  in  a  variety  ways.  The  use  of  other  coordinate  directions  may  be 
possible,  but  it  is  best  used  when  linear  transformations  of  the  random  variables  have  a  known  distributional 
form  as  is  the  case  for  normally  distributed  random  variables.  As  mentioned  above,  a  common  procedure 
is  to  partition  the  support  of  the  random  variables  and  to  apply  the  bound  on  each  of  the  partitions.  Here 
we  give  a  parametric  programming  approach  that  can  obtain  more  accurate  results  without  partitioning  the 
random  variables.  The  following  modifications  of  Algorithm  1  provide  this  basic  bound. 

Algorithm  2 

Substitute  the  following  steps  into  Algorithm  1  to  obtain 

u'($,<f>)  =  Q(lo)  +  H{4)+    £    j+<->((-)6 -(+)$). 

t>>«)Z 
Step  1' . -Solve  the  parametric  linear  program 

min{qTx  |  x(«r,0,  ar,£r)} 

for  c  €  [0,  ^Jnax  —  |r]  (or  e  €  [0,  oo)  if  £r  is  unbounded).  This  generates  a  piecewise  linear  function  /+(eer) 
with  break  points  {0,  ei, . . . ,  er},  and  with  slope  values,  qTxl+, . . . ,  qTx^. 
Then, solve  the  parametric  linear  program 

min{9rz|x(-«r,0,ar,£r)} 

for  e  G  [0,  £r  —  £™,n]  (or  €  6  [0,  oo)  if  £r  unbounded.)  We  then  obtain  a  piecewise  linear  function  /"  with 
breakpoints  {0,  e1, . . . ,  er},  and  with  slope  values,  qTx\~ ,. . . ,  qTXjT. 

Step  S!  and  Step  4':  Substitute  mint{x^  +  (i)}  for  x3+(i)  and  mint{a^~(i)}  for  x3~(i)  in  the  definitions  of 
ar+1(i)  and  a*(i)  and  substitute  maxt{xt  +  (t')}  for  x3+(i)  and  maxt{it_(i)}  for  x3~(i)  in  the  definitions  of 
/?r+1(i)  and/T(t)- 

The  changes  in  Step  l'  of  Algorithm  2  lead  to  better  bounds  if  ar  and  /9r  are  the  same  and  qT x\  <  qT  xTT. 
In  Algorithm  2,  the  approximation  obtains  as  low  a  value  as  possible  for  all  fr  for  changes  in  the  rth  direction 
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given  the  values  found  for  movement  in  previous  directions.  In  Algorithm  1,  the  approximation  just  uses  the 
extreme  values  of  £r. 

This  difference  can  be  seen  in  the  example  from  Section  4.1  which  is  illustrated  in  Figure  2.  The  dashed 
line  corresponds  to  the  function  used  in  Algorithm  1  by  using  the  extreme  values.  The  solid  line  corresponds 
to  the  functions  /+  and  /".  The  new  bound  is 

SPLU'  =  E{U' {$,$)]  =  1.449. 

We  note  that  SPLU'  is  now  below  the  EM  bound  value  of  1.625. 

The  bound  from  Algorithm  2  is  not  always  better  than  SPLU  because  the  bounds  may  change  for 
different  values  of  r,  i.e.  c*r+1  may  increase  and  /?r+1  may  decrease.  Although  this  difference  appears  to  rarely 
make  SPLU'  worse  than  SPLU  according  to  our  limited  conmputational  experience,  it  may  be  advantageous 
to  guarantee  that  a  bound  at  least  as  good  as  SPLU  is  obtained.  This  guarantee  is  accomplished  in  the 
following  modification  of  Step  1'. 

Step  1".  Solve 

min{gTx|x((^max-^)er,0,ar(^} 

T— 
=  q    x. 


Let 


a:(i)    r*(0  ifx<o, 

1 0  otherwis 


and 


I  0  otherwis< 


otherwise. 
Then  solve  the  parametric  linear  program 

min{<?Tx  |  x((€er)er,0,a^,^)} 

to  obtain  /+(eer)  as  explained  in  Step  1'. 
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This  modification  of  Algorithm  2  results  in  bounds  that  are  at  least  as  sharp  as  SPLU  and  can  still 
benefit  from  the  parametric  program  as  in  the  example  given  above.  The  key  benefit  of  the  SPLU  bound 
that  the  computational  effort  only  grows  polynomially  with  increases  in  the  number  of  random  variables  is 
maintained.  We  have  demonstrated  how  this  improvement  results  in  reduced  times  on  one  set  of  examples 
and  that  the  greatest  value  of  the  SPLU  bound  may  be  in  cases  where  the  EM  and  other  exponential  bounds 
cannot  be  reasonably  computed.  The  refinements  mentioned  above  may  allow  the  SPLU  bound  to  be  even 
more  useful  in  the  solution  of  practical  stochastic  linear  programming  problems. 
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Figure  1.  Example  with  two  random  variables. 
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Figure  2.  Parametric  linear  program  bound. 
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